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Abstract 

A  three-phase,  electrochemical  and  thermal  coupled  model  is  developed  for  valve-regulated  lead-acid  (VRLA)  batteries.  Physical 
phenomena  important  to  the  VRLA  battery  overcharge  process,  such  as  gas  generation,  transport,  and  recombination,  electrolyte 
displacement  and  capillary  flow,  and  the  venting  event  during  discharge/rest/charge,  are  incorporated  in  the  model.  The  effects  of  important 
parameters,  including  the  electrolyte  saturation  level,  interfacial  mass  transfer  coefficient  of  oxygen,  and  electrode  morphology  factor  (MF), 
are  extensively  studied.  An  overview  of  the  simulation  capabilities  of  the  present  comprehensive  model  is  provided.  ©  2002  Elsevier  Science 
B.V.  All  rights  reserved. 
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1.  Introduction 

Valve-regulated  lead-acid  (VRLA)  batteries  are  being  used 
to  power  electric  vehicles  (EV)  because  of  their  affordability 
and  reliability.  A  unique  feature  of  the  VRLA  cell  is  its 
internal  gas  cycle  in  which  oxygen  evolved  during  over¬ 
charging  of  the  positive  electrode  transports  through  a  gas- 
pore  network  to  the  negative  electrode  where  it  is  reduced 
[1,2].  As  oxygen  reactions  accompany  primary  electrode 
reactions  during  charge,  a  VRLA  battery  needs  to  be  over¬ 
charged  to  a  certain  extent  to  reach  its  full  charge  capacity.  Its 
cycle  life  can  be  greatly  affected  by  how  it  is  charged  [3]. 
While  an  advanced  charging  algorithm  for  VRLA  batteries 
may  be  achieved  through  extensive  laboratory  testing,  a 
comprehensive  first-principle-based  mathematical  model  pro¬ 
mises  to  provide  fundamental  understanding  of  the  processes 
involved  and,  thus,  serve  as  an  efficient  tool  to  search  for  an 
optimal  charging  algorithm. 

Prior  efforts  have  been  made  to  develop  mathematical 
models  for  VRLA  batteries  to  simulate  their  discharge  and 
charge  performance.  Notably,  Bernardi  and  Carpenter  [4] 
presented  an  elegant  model  of  oxygen  evolution  and  recom¬ 
bination  in  a  lead-acid  battery.  This  model,  however, 
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assumes  infinitely  rapid  transport  of  oxygen  across  the 
separator  and,  thus,  the  rate  of  oxygen  recombination  is 
limited  only  by  the  oxygen  reduction  kinetics  on  Pb  elec¬ 
trode.  As  a  result,  the  electrolyte  saturation  level  in  the 
separator  would  not  affect  the  oxygen  cycle.  In  addition, 
their  model  is  isothermal.  Newman  and  Tiedemann  [5] 
discussed  a  lumped-parameter  model  that  sought  to  clarify 
the  oxygen  evolution  and  recombination  mechanisms.  Cur¬ 
rent  and  species  concentration  distributions  were  ignored  in 
order  to  considerably  simplify  their  model.  As  a  result,  their 
model  does  not  have  resolution  to  ascertain  the  effects  of 
transport  properties  and  separator  design  on  the  internal  gas 
cycle.  Gu  et  al.  [6,7]  developed  a  comprehensive  model  that 
includes  the  following  aspects:  (1)  multi-dimensionality; 

(2)  electrochemical  kinetics;  (3)  electrolyte  flow;  (4)  species 
transport  by  diffusion,  convection  and  migration;  (5)  vari¬ 
able  porosity  of  electrodes;  (6)  variable  electrochemical  and 
transport  properties  and  (7)  cycling. 

The  objectives  of  this  paper  are  to  extend  the  lead-acid 
model  of  Gu  et  al.  [6,7]  and  with  particular  focus  on  the 
internal  oxygen  cycle  as  well  as  its  effects  on  the  full  charge 
and  overcharge  behaviors  of  VRLA  batteries.  Specifically, 
we  intend  to:  (1)  include  electrolyte  displacement  and 
capillary  flow;  (2)  account  for  oxygen  gas  transport  and 

(3)  couple  the  thermal  process  with  the  battery  model  so  that 
one  can  predict  both  electrochemical  and  thermal  behaviors 
under  non-isothermal  conditions.  While  the  present  effort  is 
directed  towards  a  single  cell,  the  developed  model  frame¬ 
work  is  expected  to  be  extendable  to  multi-cell  modules. 
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Nomenclature 

a  specific  surface  area  active  for  electrode  reac¬ 

tion  (cm2/cm3) 

c  concentration  of  species  i  in  a  phase  (mol/cm3) 

cp  specific  heat  (J/kg  K) 

ds  particle  diameter  (cm) 

D'  diffusion  coefficient  of  species  i  (cm2/s) 

EUC  electrode  utilization  coefficient  defined  as  the 
ratio  of  discharged  to  the  maximum  charge 
capacity 

F  Faraday’s  constant  (96,487  C/mol) 
h  heat  transfer  coefficient  (W/cm2  K) 

H'  Henry  constant 

iqy.ref  exchange  current  density  of  reaction  j  at  a 
reference  state  (A/cm2) 

i„j  transfer  current  density  of  reaction  j  (A/cm2) 

I  applied  current  density  (A/cm2) 

evaporation  rate  of  oxygen  at  the  electrolyte/gas 

®  3 

interface  (mol/cm  s) 

k  interfacial  mass  transfer  coefficient  at  the 

electrolyte/gas  interface  (cm/s) 
krc  relative  permeability  of  the  electrolyte  phase  in 

porous  electrode  or  separator 
K  absolute  permeability  of  a  porous  medium  (cm2) 
/  diffusion  length 

M  molecular  weight  (kg/kmol) 

p  pressure  (bar) 

Pc  valve  closing  pressure  (bar) 

PQ  valve  opening  pressure  (bar) 

q  volumetric  heat  generation  rate  (J/cm3  s) 

Q  volumetric  heat  removal  rate  from  the  cell  (J/ 
cm3  s) 

R  universal  gas  constant  (8.3143  J/mol  K) 

s  electrolyte  saturation 

S  source  term 

t  time  (s) 

t{\  transference  number  of  hydrogen  ions  with 

respect  to  the  solvent  velocity 
T  cell  temperature  (K) 

Uj  open-circuit  potential  of  electrode  reaction  j  (V) 
ve  electrolyte  velocity  vector  (cm/s) 

V  applied  voltage  (V) 

Vc  cell  volume  (cm3) 

Ve  partial  molar  volume  of  H2SO4 

Vg  gas  volume  (cm3) 

Vo  partial  molar  volume  of  H20 

Greek  symbols 
a  transfer  coefficient 

e  volume  fraction  of  a  phase,  porosity  of  a  porous 

medium 

cf)  potential  in  a  phase  (V) 

y,  S  reaction-specific  constant  in  kinetic  rate  equa¬ 
tions  in  Table  1 


rjj  surface  overpotential  of  electrode  reaction  j  (V) 

k  ionic  conductivity  of  electrolyte  (S/cm) 

kd  diffusional  conductivity  (A/cm) 

p  electrolyte  viscosity  (kg/cm  s) 

p  density  (g/cm3) 

er  conductivity  of  solid  active  material  (S/cm) 

Subscripts 
c  capillary 

e  electrolyte  phase 

eff  effective 

eg  electrolyte/gas  interface 

g  gas  phase 

s  solid  phase 

se  solid/electrolyte  interface 

sep  separator 

0  initial  value 

Superscripts 
avg  average 

eff  effective 

H  proton 

02  oxygen 


2.  Theory 

2.1.  Description  of  the  system 

As  schematically  shown  in  Fig.  1,  the  VRLA  cell  under 
consideration  consists  of  a  negative  electrode  (Pb/PbS04),  a 
positive  electrode  (Pb02/PbS04),  an  absorptive  glass  mat 
(AGM)  separator,  and  a  head  space  to  store  evolved  gas.  The 
electrolyte  is  H2S04.  Electrochemical  reactions  taking  place 
at  each  electrode  during  charge  and  discharge  are  given  as 
follows. 

Positive  electrode  (Pb02/PbS04): 

charge 

PbS04(s)  +  2H20  ?±  Pb02(s)  +  HS04~  +  3H+  +  2e~ 

discharge 

(i) 

with  oxygen  evolution 

2H20  ->  02  +  4H+  +  4e~  (2) 

and  hydrogen  recombination 

H2  2H+  +  2e-  (3) 

Negative  electrode  (Pb/PbS04): 

charge 

PbS04(s)  +  H+  +  2e~  ?±  Pb(s)  +  HS04-  (4) 

discharge 

with  oxygen  recombination 
02  +  4H+  +  4e“  2H20 


(5) 
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(-)  Pb/PbS04  Pb02  /PbS04(+) 


and  hydrogen  evolution 

2H+  +  2e  ->  H2  (6) 

Overall  cell  reaction: 

charge 

2PbS04(s)  H-  2H2O  Pb(s)  Pb02(s)  +  2H2SO4  (7) 

discharge 

The  side  reactions  (2),  (3),  (5)  and  (6)  constitute  the  internal 
oxygen  and  hydrogen  cycles  in  the  cell.  In  the  presence  of 
the  side  reactions,  a  VRLA  cell  is  a  three-phase  system 
consisting  of  the  solid  matrix,  the  liquid  electrolyte  and  a  gas 
phase.  During  charging  and  overcharging,  oxygen  is  gener¬ 
ated  at  the  Pb02/electrolyte  interface  and  may  evolve  into  the 
gas  phase  after  exceeding  its  solubility  limit  in  the  electrolyte. 
The  oxygen  can  then  be  transported,  via  the  liquid  and  gas 
phases,  from  the  positive  to  negative  electrode  where  the 
oxygen  gas  may  dissolve  back  in  the  electrolyte  and  be 
reduced  at  the  Pb/electrolyte  interface.  This  process  forms 
an  internal  oxygen  cycle  in  VRLA  cells.  Hydrogen  is  gen¬ 
erated  following  Eq.  (6)  at  the  negative  electrode  when  the 
electrode  potential  goes  too  low.  Hydrogen  recombination 
occurring  at  the  positive  electrode  as  described  by  Eq.  (3)  is 
negligible  because  of  its  poor  kinetics  [8,9].  Accumulation  of 
oxygen  and  hydrogen  in  the  gas  phase  contributes  to  the  cell 
pressure  build-up,  thereby  causing  venting. 

3.  Model  assumptions 

1 .  A  continuous  gas-pore  network  exists  in  the  cell  for  gas 
transport  and  there  is  a  head  space  right  above  the  cell 
with  a  constant  volume  to  store  excess  gas.  The  gas 
obeys  the  ideal  gas  law. 

2.  The  solid  phase  is  completely  covered  by  the  electrolyte 
film.  In  other  words,  no  mass  transfer  occurs  at  the 
active  materials/gas  interface. 


3.  Gas  convection  is  neglected  due  to  limited  gas  volume 
fraction  inside  the  cell,  leaving  the  gaseous  oxygen 
transport  by  diffusion  only. 

4.  With  the  exception  of  the  dissolved  oxygen  species, 
interfacial  chemical  and  electrical  equilibrium  exists  in 
the  electrolyte  phase  due  to  large  values  of  the  mass 
diffusivity  and  ionic  conductivity. 

5.  Hydrogen  recombination  at  the  positive  electrode  is 
neglected  due  to  its  poor  kinetics. 

4.  Mathematical  formulation 

4.1.  Governing  equations 

With  the  above  assumptions,  a  multiphase,  electrochemi¬ 
cal  and  thermal  coupled  model  for  the  system  considered  is 
derived  based  on  the  modeling  approach  presented  in 
[10,11].  Fig.  2  summarizes  the  present  model  formulation 
and  its  predictabilities.  The  battery  thermal  process  fully 
couples  with  the  electrochemical  model  via  the  heat  gen¬ 
eration  rate  and  temperature-dependent  physico-chemical 
properties.  The  derivation  of  the  momentum  equation  that 
describes  the  capillary  flow  of  the  electrolyte  is  detailed  in 
the  Appendix  A.  The  detailed  derivation  of  other  governing 
equations  is  similar  to  that  for  the  Ni-MH  battery  [12]  and, 
thus,  will  not  be  repeated  here.  The  final  model  equations  are 
summarized  in  Table  1,  along  with  source  terms  and  effec¬ 
tive  transport  properties  appearing  in  the  conservation  equa¬ 
tions  as  listed  in  Tables  2  and  3,  respectively.  Only  important 
parameters  appearing  in  Section  6  are  elaborated  in  the 
subsequent  sections. 

Electrolyte  saturation:  The  electrolyte  saturation  is 
defined  as 


£  £ 


where  s  is  the  porosity  of  an  electrode  or  separator,  £e  and 
£g  are  the  volume  fractions  of  electrolyte  and  gas  phases, 
respectively.  Oxygen  transport  from  the  positive  to  negative 
electrode  is  facilitated  by  low  saturations  of  electrolyte  (i.e. 
large  void  volume)  in  the  separator.  Culpin  and  Hayman  [13] 
showed  that  the  effective  oxygen  gas  diffusion  coefficient 
could  vary  by  three  orders  of  magnitude  if  the  saturation 
level  changes  from  90  to  70%.  Brost  [14]  further  documen¬ 
ted  the  implication  of  this  effect  to  the  field  performance  of 
VRLA  batteries.  The  electrolyte  saturation,  which  deter¬ 
mines  the  capillary  pressure,  is  a  primary  driving  force  to  the 
capillary  flow  of  the  electrolyte  [15]. 

Interfacial  mass  transfer  coefficient  of  oxygen:  The  inter¬ 
facial  evaporation  rate  of  oxygen  from  liquid  to  gas  phase  is 
determined  by 

4g  =  k(c^  -  H'cf)  (9) 

where  k  is  the  interfacial  mass  transfer  coefficient  of  the 
dissolved  oxygen  on  the  electrolyte  side  and  H'  is  the 


W.B.  Gu  et  al.  /Journal  of  Power  Sources  108  (2002)  174-184 


177 


Fig.  2.  The  thermal-electrochemical  coupled  model. 


Henry  constant.  The  difference  between  the  phase-aver¬ 
aged  concentration,  c°2,  and  the  interfacial  concentration 
of  dissolved  oxygen,  H'c °2,  provides  a  driving  force  for 
interfacial  mass  transfer  from  the  electrolyte  to  the  gas 
phase.  Using  the  concept  of  diffusion  length  proposed  by 


Wang  et  al.  [10],  the  mass  transfer  coefficient  can  further 
be  expressed  as 

D° 2 

k  =  aeg -f-  (10) 

leg 


Table  1 

Summary  of  governing  equations 


Kinetic  rate  equations  ( rjj  =  (j)s  —  (j)e  —  Uf) 


Conservation  of  charge 


Conservation  of  species 


Conservation  of  mass 


Conservation  of  momentum 


/cHyr 

Primary  reactions  (j  —  1  and  4):  inj  —  z'q/,ref  I  -q-  )  exp 

\Cref  /  L 

Oxygen  generation  and  recombination  (j  —  2  and  5): 


RT 


*lj  “exP  - 


lnj  ~  l0j,ref  [  H 
cref. 


exp 


RT 


7/  - 


.o2  V0* 


o2 


<xcjF 

expi 


,H  \  7h2 


Hydrogen  evolution  (j  =  6):  i„j  =  —  (q/.ref  (  H  I  exp 

\Cref/ 

In  electrolyte:  V  ■  (tceffV<i>e)  +  V  ■  [tcffVflnc11)]  +  =  0 

In  solid  electrode:  V  ■  (<reffV<)>s)  —  =  0 

In  electrolyte:  +  V  ■  (vecH)  =  V  ■  (£>“fVfH)  +  SH 

+  V  ■  (vec°2)  =  V  ■  (De%Vc°2)  +  S?2 
OU  c°i) 

In  gas  phase:  =  V  ■  (D°2effVc°2)  + 

dse  „  ,, 

In  electrolyte:  +  V  •  ve  =  Sj 


RT 


dt 

In  solid  electrode: 

at 


RT 


ve  =  -DcS7s  -  D£/kV or  e ^  =  V  •  Dc V5  +  V  •  Ds/K  V  ^  +  sS sv ) 

d (P£pT)  . 


dr 


=  q~Q 


Conservation  of  thermal  energy 
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Table  2 

Source  terms  in  conservation  equations 


Si 

si 


Terms 

Positive  electrode 

Separator 

Negative  electrode 

5 

at  ini  +  U2i„2 

0 

^4  ln4  “1“  #5  ln5  “1“  ^6  ln6 

SH 

3-2/“  1  -  /“ 

1  —  2ft  1  - 

_  a\in\  +  _  aiini 

2  F  F 

0 

a4in  4+  _  (a5in5  +  d6in6) 

zr  r 

_/02 

eg 

—  fl5i„5  -  2°2 

J°  2 
eg 

D°2 

k(c°2  -H'c°2),k  =  aes-^, 

leg 

[Ve(3  -  2r° )  -  2V0]  +  [Ve(2  -  2t°+)  -  Vo] 

( M?b02  M] pbS04 \  alln\ 


V  Ppb02  PPbS04  7  2Z7 


Ve(l  -  2t°+)  ^  +  [Ve(2  -  2t°+ )  -  Vo]  ^  +  Ve(2  -  2f°  )  ^ 


/  ^PbS04 
V  PpbSC>4 


Mpb\ 

Ppb  /  2F 


dV  -  /  [*e"V4>e  ■  V0e  +  ■  V0JdV 

V c  ./V, 


where  «eg  is  the  specific  area  of  the  electrolyte/gas  interface, 
D the  diffusion  coefficient  of  the  dissolved  oxygen  in  the 
electrolyte,  and  /eg  is  the  microscopic  diffusion  length  and  its 
value  is  proportional  to  the  electrolyte  film  thickness  [10]. 
Oxygen  recombination  at  the  negative  electrode  is  limited 
by  the  interfacial  mass  transfer  coefficient.  A  thin  film 
condition  in  the  negative  electrode  will  facilitate  the  oxygen 
transport  through  the  electrolyte/gas  interface  and.  hence, 
the  oxygen  recombination. 

Electrode  morphology  factor  ( MF ):  A  principal  geometric 
parameter  for  characterizing  electrode  performance  is  the 
specific  area  active  for  electrode  reaction  j,  which  is,  in  turn, 
dependent  upon  the  porosity  and  particle  size  as  well  as  the 
local  state  of  charge  (SOC),  namely, 

asj  =  aSj(l  —  EUC’),  for  discharge  i  =  Pb02  orPb 
«s;  =  as/EUO ,  for  charge  i  =  Pb02  or  Pb 


Table  3 

Effective  transport  properties  in  conservation  equations 


Symbol 

Expression 

Reference 

Ele5K 

Bruggeman  relation 

eff 

kd 

RT Keff ,  „ 

P  (2<+  —  1) 

[6] 

Bruggeman  relation 

nH 

ueff 

e‘5Dh 

Bruggeman  relation 

n°2 

ue,elf 

Bruggeman  relation 

n° 2 

^g,eff 

Dc 

f(s)D°2 

/ey/2  Kkre 

-g{k) 

empirical  correlation 

Eq.  (A.  10) 

DfJK 

j/  \  Kkre 

— (jJ(s) - 

ft 

180(1  -  s)2 

Eq.  (A.  11) 

K 

Kozeny-Carman  equation 

for  primary  reactions  and 

asj  =  fls,.(l  -  EUC’) ,  i  =  Pb02  or  Pb  (12) 

for  oxygen  and  hydrogen  reactions.  Here,  as  denotes  the 
specific  area,  with  subscript  i  denoting  the  active  materials 
Pb02  at  the  positive  electrode  and  Pb  at  the  negative 
electrode,  respectively  and  EUC  represents  the  electrode 
utilization  coefficient  defined  as  the  ratio  of  used  charge 
capacity  to  the  maximum  charge  capacity  at  fully-charged 
state,  Eqs.  (11)  and  (12)  describe  the  amount  of  active  area 
during  discharge  or  charge,  and  the  exponent  £,  is  a  MF 
describing  the  rate  of  change  in  the  specific  area. 

Cell  pressure  and  gas  partial  pressures :  The  average  cell 
pressure  can  be  calculated  by  summing  up  partial  pressures 
of  all  gases  present  in  the  gas  phase: 

P7  =  Yy ’  *  =  °2  and  N2  (13) 

i 

where  the  nitrogen  gas  partial  pressure  is  determined  from 
the  initial  condition  and  remains  unchanged  until  venting 
takes  place.  The  oxygen  and  hydrogen  partial  pressures  are 
calculated  using  the  ideal  gas  law: 

p‘ =  cfsRT,  i  =  02  and  H2  (14) 

The  average  oxygen  and  hydrogen  concentrations  in  the  gas 
phase  are  evaluated  by 

Cvg  7^  /  dV,  i  =  02  and  H2  (15) 

Pc  Jvc 

where  Vc  is  the  entire  cell  volume  including  that  of  the  head 
space. 

In  a  VRLA  battery,  the  pressure  is  regulated  through  a 
venting  valve.  Whenever  the  internal  gas  pressure  becomes 
greater  than  the  valve  preset  opening  pressure,  PQ ,  gas  is 
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released  into  the  ambient  until  the  internal  pressure 
drops  below  the  valve  preset  closing  pressure,  Pc.  Assum¬ 
ing  that  the  time  of  gas  relief  at  pressures  higher  than 
the  opening  pressure  of  the  valve  is  negligible,  a  single 
venting  event  can  be  simply  modeled  by  the  following 
procedure: 

if  PgVg  <  P0,  thenftg  =  ng  and  p™g  =  p'd^g 

n  P.  (16) 

if  p™g  >  P0,  then  n„  =  gavg  and/^vg  =  Pc 
Pg 

where  n„  is  the  number  of  moles  in  gas.  All  gas  partial 
pressures  and  their  local  concentrations  decrease  accord¬ 
ingly  with  venting. 

4.2.  Initial  and  boundary  conditions 


That  is,  there  is  no  electrolyte  flow  initially.  Electrolyte  starts 
to  flow  after  electrode  reactions  take  place,  resulting  in 
changes  in  volume  fractions  of  solid  and  electrolyte,  and 
thus,  changes  in  the  electrolyte  saturation  level,  porosity,  and 
permeability. 

Eq.  (23)  is  used  to  evaluate  the  initial  electrode  perme¬ 
ability  given  the  permeability  of  separator.  At  electrode/ 
separator  interfaces,  there  are  sharp  changes  of  electro¬ 
lyte  saturation  and  porosity.  Eq.  (23)  shows  that  the  velo¬ 
city  induced  by  the  saturation  gradient  is  balanced  by 
the  velocity  caused  by  the  gradient  of  ( e/K)v 2  at  equili¬ 
brium. 


5.  Numerical  procedures 


To  close  the  mathematical  system,  initial  and  boundary 
conditions  for  the  primary  variable  are  necessary.  The 
potentials  in  the  solid  and  electrolyte  phases  are  governed 
by  Poisson  equation  without  time  derivatives,  hence,  their 
initial  conditions  are  not  necessary.  The  appropriate  bound¬ 
ary  conditions  for  <fis  and  cj)e  are 

<^s  =  Oor  V  (17) 

on  the  current  collectors  for  a  given  voltage,  or 


To  numerically  solve  the  foregoing  model  equations,  the 
general  computational  fluid  dynamics  (CFD)  methodology 
is  adopted  in  the  present  work.  The  key  to  a  successful 
application  of  CFD  techniques  is  that  all  the  equations  in  the 
present  battery  model  can  be  cast  into  a  general  form  of 
convection-diffusion  type,  i.e. 

+  V  •  (v<2>)  =  V  •  (r#V$)  +  S*  (24) 

9t  .  convention  diffusion  source 

accunulation 


S  J 
=  1 
on 


(18) 


for  a  prescribed  current  density.  Everywhere  else  on  the  cell 
boundaries,  one  has 


dn  dn 


=  0 


(19) 


where  n  is  the  unit  vector  normal  to  the  boundary  surface. 
The  concentration  initial/boundary  conditions  are  given  by 


dc' 

c‘  =  c‘0,  at  r  =  0;  g^  =  0 


(20) 


The  temperature  initial/boundary  conditions  can  be  written 
as 


T=T0,  att  =  0:  -Ac—  =  h(T-Ta)  (21) 

where  /.c  is  the  thermal  conductivity  of  the  case  material,  and 
h  the  overall  convective  heat  coefficient  that  may  include  the 
effects  of  other  heat  transfer  modes  and  7a  is  the  ambient 
temperature. 

The  initial  conditions  for  electrolyte  velocity  and  satura¬ 
tion  level  is  determined  by  assuming  that  the  system  is 
initially  at  equilibrium,  i.e. 


where  $  is  the  general  variable  to  be  solved,  rcI,  the  diffusion 
coefficient,  and  S  is  the  source  term  which  includes  all  terms 
that  cannot  be  included  in  the  previous  terms. 

The  general  differential  equation  is  then  discretized  by  the 
control  volume-based  finite  difference  method  of  Patankar 
[16]  and  the  resulting  set  of  algebraic  equations  is  iteratively 
solved.  The  numerical  solver  for  the  general  differential 
equation  can  be  repeatedly  applied  for  each  scalar  variable 
over  a  control  volume  mesh.  The  solution  of  each  of  the 
components  of  the  electrolyte  momentum  equation  is, 
however,  obtained  in  a  staggered  control  volume  using 
semi-implicit  method  for  pressure-linked  equations  revised 
(SIMPLER)  algorithm  (see  [16]  for  details).  The  rectangular 
physical  domain  is  divided  by  either  a  uniform  or  non- 
uniform  grid  consisting  of  m  horizontal  and  n  vertical  grid 
lines.  Stringent  numerical  tests  are  performed  to  ensure  that 
the  solutions  are  independent  of  the  grid  size  and  time 
step.  The  equations  are  solved  as  a  simultaneous  set,  and 
convergence  is  considered  to  be  reached  when  the  relative 
error  in  each  field  between  two  consecutive  iterations  is  less 
than  10“5. 


6.  Results  and  discussion 


ve  =  0,  atf  =  0 
and 

/  e  \  t/2 

DcS7s  +  De/kX7  J  =0,  atf  =  0 


(22)  Fig.  3  compares  the  model  predictions  with  the  experi¬ 

mental  data  from  charging  a  VRLA  cell  at  C/10  with  40% 
overcharge  [4],  showing  reasonable  agreement  in  predicting 
both  the  cell  potential  and  pressure.  Specifically,  the  model 
successfully  captures  two  key  features  observed  in  charging 
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Fig.  3.  Comparisons  with  experimental  data. 


of  a  VRLA  cell — the  potential  rollover  and  pressure  stabi¬ 
lization  as  overcharge  proceeds. 

In  the  following  figures,  the  unique  capabilities  of  the 
present  model,  namely,  (1)  simulating  oxygen  gas  transport; 
(2)  venting  event  and  (3)  capillary  flow,  will  be  illustrated. 
Numerical  simulations  are  performed  to  study  the  effects  of 
important  parameters  such  as  the  electrolyte  saturation  level, 
interfacial  mass  transfer  coefficient  of  oxygen,  and  MF  of  an 
electrode.  The  effects  of  venting  event,  charging  algorithm, 
and  capillary  flow  are  also  explored.  The  simulated  battery 
has  a  nominal  capacity  of  70  Ah.  It  is  configured  with  four 
cells  in  series,  so  its  open-circuit  voltage  is  approximately 
8.6  Vat  100%  SOC.  In  the  following  studies,  the  response  of 


Time  (h) 

Fig.  4.  Effect  of  electrolyte  saturation  in  the  separator. 


the  battery  charged  at  a  constant  current  of  25  A  (C/2.8)  will 
be  explored. 

Fig.  4  shows  the  effect  of  electrolyte  saturation  level  on 
the  cell  performance.  Three  values  of  acid  saturation  level  in 
the  separator  are  considered:  83,  88  and  93%,  while  85% 
acid  saturation  is  assumed  for  the  positive  and  negative 
electrodes.  Fig.  4  shows  that  the  acid  saturation  level  does 
not  have  significant  effect  on  the  battery  performance  until 
oxygen  reactions  come  into  play  when  approaching  over¬ 
charge.  A  smaller  value  of  acid  saturation  level  (i.e.  a  larger 
void  volume  for  the  gas  phase  and,  hence,  larger  effective 
oxygen  gas  diffusion  coefficient)  yields  a  more  pronounced 
potential  rollover  and  a  lower  stabilized  pressure,  indicating 
that  oxygen  transport  from  the  positive  electrode,  where  it  is 
generated,  to  the  negative  electrode,  where  it  is  recombined, 
through  the  separator  indeed  plays  an  important  role  in  the 
operation  of  VRLA  batteries.  The  increase  in  oxygen  recom¬ 
bination  rate  at  lower  acid  saturation  has  been  seen  to  be 
important  in  designing  an  optimal  charging  algorithm  [3]. 

The  sharp  increase  in  the  cell  potential  at  the  end  of  charge 
occurs  when  the  negative  electrode  is  completely  charged 
and  available  oxygen  is  not  sufficient  to  provide  a  balanced 
amount  of  current  to  that  of  the  positive  electrode.  There¬ 
fore,  it  is  an  indication  of  oxygen  mass  transport  limitation, 
as  reflected  by  the  sequential  occurrence  of  the  potential 
jumps  in  Fig.  4.  This  can  be  further  illustrated  by  Fig.  5, 
which  shows  the  effect  of  oxygen  transport  from  the  positive 
to  the  negative  electrode.  In  Fig.  5,  the  solid  line  represents 
the  base  case  with  finite  values  of  the  interfacial  mass 
transfer  coefficient  k  and  oxygen  gas  diffusion  coefficient 
D° 1 ,  and  the  dashed  line  stands  for  the  case  with  infinitely 
large  values  of  k  and  //,’-  (i.e.  no  mass  transfer  resistance  is 
present  for  oxygen).  As  a  result,  all  oxygen  generated  in  the 
positive  electrode  is  consumed  in  the  negative  electrode.  No 
cell  pressure  build-up  due  to  the  accumulation  of  oxygen  is 
seen.  The  reaction  currents  are  always  balanced  so  that  the 
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Fig.  5.  Effect  of  oxygen  gas  transport. 
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cell  potential  gradually  approaches  a  constant  and  no  other 
side  reactions  come  into  play.  Fig.  5  only  gives  a  combined 
effect  of  k  and  D°2.  Simulations  performed  to  separate  the 
two  effects  suggest  that  while  the  oxygen  transport  in  the  gas 
phase  is  greatly  affected  by  the  acid  saturation  level  through 
the  effective  oxygen  gas  diffusion  coefficient,  the  interfacial 
mass  transfer  of  the  dissolved  oxygen  is  more  crucial  to  the 
oxygen  recombination  in  the  negative  electrode. 

Fig.  6  shows  the  effect  of  MF  that  accounts  for  the 
passivation  effect  of  PbSO.4  on  the  specific  electroactive 
area,  as  shown  by  Eqs.  (11)  and  (12).  During  charging,  a 
smaller  value  of  MF  yields  a  slower  change  in  the  specific 
area.  Consequently,  the  surface  overpotential  that  drives  the 
electrode  reaction  is  smaller,  and  the  potential  of  the  positive 
electrode  remains  comparatively  low  and,  thus,  oxygen 
evolution  takes  place  later.  On  the  other  hand,  a  larger  value 
of  MF  yields  a  higher  surface  overpotential  and  the  oxygen 
reactions  occur  earlier  and  become  dominant.  As  a  result,  the 
final  jump  in  the  cell  potential  disappears  as  the  negative 
electrode  has  not  yet  achieved  fully-charged  state. 

Fig.  7  compares  the  cell  performance  between  cases  with 
and  without  venting.  The  valve-opening  and  closing  pres¬ 
sures  are  0.17  and  0.10  atm  (gauge  pressure),  respectively.  It 
is  shown  that  with  venting  the  cell  pressure  is  windowed 
within  the  preset  pressures.  The  cell  potential  increases 
faster  with  venting  due  to  the  faster  conversion  of  PbS04 
to  Pb  at  the  negative  electrode,  as  the  primary  reaction 
current  remains  higher  when  oxygen  generated  in  the  posi¬ 
tive  electrode  releases  to  the  ambient  and  not  recombined  in 
the  negative  electrode  to  consume  a  portion  of  total  current. 

Fig.  8  shows  the  effect  of  hydrogen  evolution  in  the 
negative  electrode  on  the  evolutions  of  cell  potential  and 
pressure.  Since  hydrogen  recombination  in  the  positive 
electrode  is  kinetically  hindered,  hydrogen  gas,  once  gen¬ 
erated,  immediately  contributes  to  the  cell  pressure  build¬ 
up.  The  cell  potential  responses  to  the  hydrogen  evolution 


with  a  large  step  rise  followed  by  gradual  decrease  due  to 
oxygen  recombination  as  overcharge  proceeds.  Fig.  8  also 
shows  the  effect  of  venting  in  the  presence  of  hydrogen 
evolution.  The  difference  in  the  cell  potential  between  cases 
with  and  without  venting  is  substantial.  In  particular,  no 
gradual  decrease  in  the  cell  potential  is  seen  because  sig¬ 
nificant  amount  of  oxygen  generated  at  the  positive  elec¬ 
trode  is  vented  rather  than  recombined  at  the  negative 
electrode.  In  other  words,  oxygen  recombination  is  ineffi¬ 
cient,  causing  hydrogen  excursions  to  remain  high.  Such 
behavior  of  cell  potential  is  unique  for  the  VRLA  battery 
experiencing  significant  gassing  and  venting  [17]. 

Fig.  9  displays  the  temperature  evolutions  of  the  battery 
under  various  charging  algorithms.  Three  charging  algo¬ 
rithms  are  considered:  continuous  current  (CC),  current 
interrupt  (Cl),  and  control  charge.  The  CC  charging  algo¬ 
rithm  involves  an  initial  high  rate,  followed  by  a  number  of 
step-wise  decreased  charging  rates,  and  finally  a  fixed  low 
rate  that  may  be  higher  than  the  preceding  one.  Under  Cl 


Fig.  8.  Effect  of  hydrogen  evolution. 
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Fig.  11.  Effect  of  capillary  flow  on  the  distributions  of  electrolyte 
saturation  during  C/1  discharge. 


charging  algorithm  the  battery  is  charged  at  a  nearly  constant 
high  rate  and  then  pulse-charged  at  a  relatively  low  rate  near 
the  fully-charged  state  [3].  The  control  charge  is  similar  to 
the  Cl  but  more  complex. 

The  battery  is  exposed  to  the  ambient  and  is  cooled 
primarily  by  free  air  convection.  Significant  differences  in 
the  battery  temperature  occur  during  overcharge.  Essen¬ 
tially,  CC  yields  a  continuous  increase  in  temperature  with 
possibility  of  thermal  runaway  at  the  end  of  charging.  Cl 
caps  the  temperature  below  37  °C  and  not  only  achieves  full 
capacity  but  also  ensure  safe  operation  of  the  battery.  The 
temperature  evolution  in  the  control  charge  process  also 
remains  low  with  longer  charging  time.  The  battery  experi¬ 
ences  a  decrease  in  its  temperature  in  all  cases  in  accordance 
with  the  high  to  low  rate  transition  involved  in  a  charging 
algorithm.  Fig.  9  demonstrates  a  unique  capability  of  the 
present  thermal-electrochemical  coupled  model. 


Fig.  10  shows  the  effect  of  capillary  flow  on  the  battery 
voltage  and  pressure  during  a  simple  discharge/rest/charging 
cycle.  With  a  preset  cut-off  voltage,  the  simulation  con¬ 
sidering  capillary  flow  yields  about  10%  more  in  the  dis¬ 
charge  capacity,  due  to  the  fact  that  capillary  flow  enhances 
the  electrolyte  supply  from  the  separator  to  the  positive 
electrode  where  the  electrolyte  tends  to  be  depleted  during 
discharge.  The  difference  in  cell  potential  between  cases 
considering  and  ignoring  capillary  flow  is  negligible  during 
charging.  The  difference  in  the  cell  pressure  is  also  small 
because  most  changes  are  masked  by  the  venting  events  that 
occur.  However,  substantial  differences  in  the  distributions 
of  acid  saturation  level  are  observed  from  Fig.  1 1  for 
discharge  and  Fig.  12  for  charging.  While  the  case  con¬ 
sidering  capillary  flow  predicts  a  more  even  distribution  of 
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acid  saturation  level  which  reflects  the  reality,  the  case 
ignoring  capillary  flow  yields  highly  non-uniform  profiles 
of  acid  saturation  level.  These  predictions  illustrate  the 
importance  of  including  capillary  flow  in  battery  models. 
The  effect  is  expected  to  be  significant  in  predicting  battery 
cycle  life  where  there  are  a  large  number  of  deep  charge/ 
discharge  cycles. 


7.  Conclusion 

This  paper  provides  an  overview  of  the  simulation  cap¬ 
abilities  of  the  present  multiphase,  electrochemical  and 
thermal  coupled  model  for  VRLA  batteries.  Numerical 
simulations  were  carried  out  to  show  the  effects  of  the 
electrolyte  saturation  level,  interfacial  mass  transfer  coeffi¬ 
cient  of  oxygen,  MF,  venting,  charging  algorithm,  and 
capillary  flow  on  the  battery  performance.  Accurate  estima¬ 
tion  of  the  parameters  that  govern  these  phenomena  would 
allow  this  comprehensive  model  to  be  used  in  designing 
proper  charging  algorithms  for  better  performance  and  cycle 
life. 
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Appendix  A.  Derivation  of  the  momentum  equation 
for  capillary  flow 

Conservation  of  mass  over  the  liquid  phase  yields 

^  +  V-ve  =  Sev  (A.1) 

where  the  term  Sve  represents  electrolyte  displacement  dur¬ 
ing  battery  operation  resulting  from  the  change  in  partial 
molar  volumes  between  the  reactants  and  products.  The  flow 
of  electrolyte  in  porous  media  is  usually  described  by 
Darcy’s  law,  i.e, 

ve  =  -  V pe  (A.2) 

P 

where  pe  is  the  pressure  in  the  electrolyte,  K  the  absolute 
permeability  of  the  porous  medium  like  electrode  and 
separator,  kle  the  relative  permeability  of  the  phase  of 
electrolyte,  and  //  is  the  dynamic  viscosity  of  the  electrolyte. 
Since  the  electrolyte  pressure  differs  from  the  gas  pressure 
by  a  capillary  pressure  [15],  namely, 


(A.3) 


Eq.  (A.2)  can  be  rewritten  as 


ve 


Kkrc 

p 


Vpc 


Kkre 

P 


Vpg 


(A.4) 


The  first  term  on  the  right-hand  side  of  Eq.  (A.4)  describes 
the  electrolyte  flow  driven  by  capillary  forces,  whereas  the 
second  term  represents  the  electrolyte  displacement  by  the 
gradient  in  the  gas  pressure.  The  first  term  is  dominant  when 
the  thickness  of  an  unsaturated  cell  is  small  and  as  a  result, 
the  gas  phase  pressure  is  nearly  uniform.  Neglecting  the 
gradient  in  gas  phase  pressure,  one  has 


Ve 


Kkrc 

P 


Vpc 


(A.5) 


The  absolute  permeability  of  a  porous  medium  can  be 
evaluated  by  the  Kozeny-Carman  relation 


180(1  —  e)2 


(A.6) 


where  s  is  the  porosity  of  the  porous  medium  and  ds  is  the 
diameter  of  spherical  particles  making  up  it. 

The  relative  permeability  of  the  electrolyte  phase  is  a 
function  of  electrolyte  saturation  [15] 

krc  =  s"  (A.7) 


where  the  index  n  is  an  empirical  value  varying  between  1  and  3. 

The  two-phase  capillary  pressure  can  be  expressed  in 
terms  of  a  Leverett  function  J(s)  [15],  namely, 


Making  use  of  Eqs.  (A.6-A.8),  and  the  definition  of  elec¬ 
trolyte  saturation  Eq.  (8),  Eq.  (A.5)  becomes 


ve  =  ~Dc\7s  -  De/k\7 


© 


1/2 


(A.9) 


where 


(A.  10) 


and 

n  T(^Kk « 

De/K  =  -aj[s) - 


(A.l  1) 


Dc  represents  capillary  diffusion  coefficient  and  has  a  unit  of 
cnr/s.  Substituting  Eq.  (A.9)  into  Eq.  (A.l),  one  has 


=  V  •  DcVs  +  V  •  De/KV (|)  1/2+(5ey  +  <) 

(A.  12) 
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